model.name = "FXaExample1" # root name of modeling files
workDir = "~/bugsParallel/FXaExample1" # change to match the path to the directory containing FXaExample1.R and FXaExample1.txt
toolsDir = "~/bugsParallel/bugsTools" # change to match your setup
bugsDir = paste(Sys.getenv("HOME"),"/Program Files/WinBUGS14/",sep="")
nSlaves = 4

library(Rmpi)
library(rlecuyer)
library(coda)
library(lattice)
source(paste(toolsDir,"/bgillespie.utilities.R",sep="")) # These are a few of my favorite things
source(paste(toolsDir,"/bugsParallel-1.0.R",sep=""))
source(paste(toolsDir,"/bugs.tools.R",sep="")) # A few BUGS-specific utilities

set.seed("10271998") # not required but forces consistent results
setwd(workDir)
###################################################################
# Data management

xdata = read.csv("fxa.data.avg.csv")

# create WinBUGS data set
bugsdata = list(
	nobs = nrow(xdata),
	cobs = xdata$cavg,
	FXa = xdata$fxa.inh.avg)

# create initial estimates
bugsinit = function() list(
	Emax = runif(1,40,100),
	logEC50 = rnorm(1,log(100),0.4),
	gamma = 10*rbeta(1,0.25,5),
	sigma = exp(rnorm(1,log(5),0.2)))

# specify what variables to monitor
parameters = c("Emax","EC50","gamma","sigma","FXaPred")

# specify the variables for which you want history and density plots
parameters.to.plot = c("deviance","Emax","EC50","gamma","sigma")

################################################################################################
# run WinBUGS

n.chains = nSlaves
n.iter = 60000
n.burnin = 10000
n.thin = 50

mpi.spawn.Rslaves(nslaves = nSlaves) # launches multiple R slave processes on available processors 
mpi.setup.rngstream() # initializes parallel random number generation

if(.Platform$OS.type == "windows"){
	bugs.fit <- bugsParallel(data=bugsdata,inits=bugsinit,
		parameters.to.save=parameters,model.file=paste(getwd(),"/",model.name,".txt",sep=""),
		n.chains=n.chains,n.iter=n.iter,n.thin=n.thin,n.burnin=n.burnin,refresh=10,clearWD=T,
		bugs.directory = bugsDir)
}else if(.Platform$OS.type == "unix"){
#	wineBin = "/Applications/Darwine/Wine.bundle/Contents/bin"
	wineBin = "/usr/local/MacPorts/bin"
	bugs.fit <- bugsParallel(data=bugsdata,inits=bugsinit,
		parameters.to.save=parameters,model.file=paste(getwd(),"/",model.name,".txt",sep=""),
		n.chains=n.chains,n.iter=n.iter,n.thin=n.thin,n.burnin=n.burnin,refresh=10,clearWD=T,
		bugs.directory = bugsDir, useWINE=T, WINE=paste(wineBin,"/wine",sep=""),
		newWINE=T, WINEPATH=paste(wineBin,"/winepath",sep=""))
}
 
# save scripts, data and results to a directory
save.model(bugs.fit,model.name)

# convert MCMC results to formats suitable for post-processing
sims.array = aperm(array(unlist(bugs.fit),dim=c(nrow(bugs.fit[[1]]),ncol(bugs.fit[[1]]),length(bugs.fit)),
	dimnames=c(dimnames(bugs.fit[[1]]),list(NULL))),c(1,3,2))
posterior = array(as.vector(sims.array),dim=c(prod(dim(sims.array)[1:2]),dim(sims.array)[3]),
	dimnames=list(NULL,dimnames(sims.array)[[3]]))

################################################################################################
# posterior distributions of parameters

# open graphics device
pdf(file = paste(model.name,"/",model.name,".plots.pdf",sep=""),width=6,height=6)

# create history, density and Gelman-Rubin-Brooks plots, and a table of summary stats
ptable = parameter.plot.table(sims.array[,,unlist(sapply(c(paste("^",parameters.to.plot,"$",sep=""),
	paste("^",parameters.to.plot,"\\[",sep="")),grep,x=dimnames(sims.array)[[3]]))])

write.csv(signif(ptable,3),paste(model.name,"/",model.name,".summary.csv",sep=""))

################################################################################################
# posterior predictive distributions

pred = posterior[,grep("FXaPred\\[",dimnames(posterior)[[2]])]
x1 = xdata
x1$type =rep("observed",nrow(x1))
x2 = rbind(x1,x1,x1)
x2$fxa.inh.avg = as.vector(t(apply(pred,2,quantile,probs=c(0.95,0.5,0.05))))
x2$type = rep(c("5%ile","median","95%ile"),ea=nrow(x1))
x1 = rbind(x1,x2)
x1 = x1[order(x1$type,x1$cavg),]

xyplot(fxa.inh.avg~cavg,x1,groups=type,panel=panel.superpose,pch=c(NA,NA,NA,1),
	type=c("l","l","l","p"),lty=c(3,3,1,0),col=c("red","red","blue","black"),lwd=3,
	scales=list(cex=1),par.strip.text=list(cex=1),
	strip = function(...) strip.default(..., style = 1),
	xlab=list(label="time-averaged plasma drug concentration",cex=1.2),
	ylab=list(label="time-averaged FXa inhibition (%)",cex=1.2))

dev.off() # close graphics device

mpi.close.Rslaves() # stop R slave processes

